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Abstract.  Motivated  by  recent  efforts  to  mitigate  blast  loading  using  energy-absorbing  materials,  this 
paper  investigates  the  uniaxial  crushing  of  cellular  media  in  sandwich  construction  under  impulsive 
pressure  loading.  The  cellular  core  is  modeled  using  a  rigid,  perfectly-plastic,  locking  idealization,  as 
in  previous  studies,  and  the  front  and  back  faces  are  modeled  as  rigid,  with  pressure  loading  applied  to 
the  front  face  and  the  back  face  unrestrained.  Predictions  of  this  analytical  model  show  excellent 
agreement  with  explicit  finite  element  computations,  and  the  model  is  used  to  investigate  the  influence 
of  the  mass  distribution  between  the  core  and  the  faces.  Increasing  the  mass  fraction  in  the  front  face  is 
found  to  increase  the  impulse  required  for  complete  crushing  of  the  cellular  core  but  also  to  produce 
undesirable  increases  in  back-face  accelerations.  Optimal  mass  distributions  are  investigated  by 
maximizing  the  impulse  capacity  while  limiting  the  back-face  accelerations  to  a  specified  level. 

Keywords:  Blast  mitigation,  aluminum  foam,  shock  wave,  finite  element  analysis. 

PACS:  46.40.Cd,  62.50.+p,  83.60.Uv. 


BACKGROUND 

Cellular  materials  such  as  metal  foams  and 
honeycombs  are  being  considered  in  a  wide  variety 
of  structural  applications  because  of  their  capacity 
to  absorb  impact  energy.  Surprisingly,  however, 
their  use  under  blast  loading  has  often  led  to 
enhancement,  rather  than  mitigation,  of  blast 
effects.  Experiments  by  Hanssen  et  al.  [1]  showed 
that  increased  upswing  results  from  the  addition  of 
an  aluminum  foam  layer  to  the  face  of  a  massive 
“pendulum”  subjected  to  blast  loading.  Nesterenko 
[2]  noted  that  in  these  experiments,  the  blast 
impulse  is  imparted  primarily  to  a  lightweight  plate 
covering  the  foam  layer,  leading  to  significantly 
higher  kinetic  energy  than  if  the  same  impulse 
were  imparted  directly  to  the  more  massive 
pendulum.  Xue  and  Hutchinson  [3]  noted  a  similar 


effect  in  a  computational  study  of  blast  loading  on 
sandwich  plates,  in  which  the  kinetic  energy 
imparted  to  a  sandwich  plate  was  observed  to  be 
greater  than  for  a  solid  plate  of  the  same  mass.  In 
spite  of  this,  it  was  found  that  deflections  of 
sandwich  plates  could  be  significantly  less  than  for 
the  corresponding  solid  plate.  Xue  and  Hutchinson 
considered  front  and  back  face  sheets  with  equal 
mass  but  suggested  that  further  reductions  in 
deflections  might  be  achieved  by  increasing  the 
mass  fraction  in  the  face  sheet  near  the  blast. 

ANALYTICAL  MODEL 

Motivated  by  these  observations,  an  analytical 
model  is  developed  in  this  paper  to  investigate  the 
influence  of  mass  distribution  on  the  uniaxial 
crushing  of  cellular  material  sandwiched  between 
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rigid  layers.  The  cellular  core  material  is 
represented  by  the  simplified  stress-strain 
relationship  shown  in  Fig.  1(b),  originally  proposed 
by  Reid  and  Peng  [4]  for  modeling  crushing  of 
wood  and  subsequently  applied  to  cellular  metals 
in  a  number  of  studies  (e.g.,  [1,5]).  Arbitrary 
masses  of  the  front  and  back  faces  are  permitted, 
and  a  pressure  pulse  p{t)  is  applied  to  the  front  face 
with  the  back  face  unrestrained.  This  sandwich 
model  is  a  generalization  of  that  in  [1],  which 
considered  a  fixed  back  face,  and  of  that  in  [5], 
which  considered  front  and  back  faces  of  equal 
mass  with  blast  loading  represented  by  an  initial 
velocity  imparted  to  the  front  face. 

A  strip  of  sandwich  panel  with  unit  cross- 
sectional  area  is  considered,  with  total  mass  given 
by  w  =  +^2 ,  where  and  are  the 

uncompressed  density  and  thickness  of  the  cellular 
core,  and  and  are  the  areal  densities  of  the 
front  and  back  faces.  The  acceleration  of  the  center 
of  mass,  denoted  ,  follows  directly  from 
application  of  Newton’s  second  law  to  the  strip: 

p{t)  =  mUg  (1) 

Provided  the  applied  pressure  is  sufficiently  high, 
densification  of  the  cellular  core  commences  at  the 
front  face,  and  a  densification  front  propagates 
through  the  core.  By  conservation  of  mass,  the 
density  of  the  compressed  core  material  is 
According  to  the  simplified  model  of 
Fig.  1(b),  the  compressed  core  material  moves  as  a 
rigid  body  with  the  same  velocity  as  the  front  face, 
denoted  ,  while  the  uncompressed  core  material 
moves  as  a  rigid  body  with  the  velocity  of  the  back 
face,  t/2  •  The  stress  just  ahead  of  the  densification 
front  is  cTq,  and  application  of  Newton’s  second 


law  to  the  material  ahead  of  the  densification  front 
then  yields  the  following  equation: 

O-Q  =(Po^  +  '”2)«2  (2) 

where  x  denotes  the  thickness  of  the  uncompressed 
core  material,  and  the  thickness  of  the  densification 
front  itself  is  assumed  to  be  negligible.  By  forming 
and  differentiating  an  expression  for  x^. ,  the 
distance  of  the  center  of  mass  from  the  back  face,  it 
follows  that 

Xq  ={sjm)^^m^+ p^{l^-x)\x-p^x'"]  (3) 

Eqs.  (1)  -  (3)  can  then  be  combined  through  the 
relation  +  x^.  to  yield  the  following 

nonlinear  ordinary  differential  equation  for  x: 

-E^[m,+ p^{i^-x)]x  +  E^p^x^ 

=  p{t)  -G^m/  (PqX  +  ^2  ) 

Eq.  (4)  can  be  integrated  numerically  with  initial 
conditions  x(0)  =  f  ^  and  x(0)  =  0  .  A  triangular 
pressure  pulse  is  considered,  as  shown  in  Fig.  1(c), 
with  total  impulse  denoted  4.  The  following 

symbols  are  introduced  to  denote  the 

nondimensional  peak  pressure  and  total  impulse: 


The  following  symbols  denote  the  mass  fractions  in 
the  core  and  in  the  front  and  back  faces: 

t/q  =  yOpf  Q  !  m  \  m  (6) 


Figure  1.  Analytical  model  definition:  (a)  Strip  of  sandwich  panel  with  partially  compacted  core;  (b)  Stress  vs. 
volumetric  strain  relationship  for  core  material;  (c)  Triangular  pressure  pulse  applied  to  front  face. 
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Table  1.  Parameters  of  eomputational  simulations. 


Case 

% 

72 

Po 

4 

pendulum 

0.0125 

0.0125 

0.975 

10 

0.015 

sandwich 

0.25 

0.5 

0.25 

10 

1 

COMPARISON  WITH  COMPUTATIONS 


defined  as  v^-u^l  and  V2  where 

^00  -  4  /  is  velocity  of  the  center  of 

mass.  Due  to  the  small  mass  of  the  front  face, 
much  larger  nondimensional  front-face  velocities 
are  observed  in  the  “pendulum”  case,  despite  the 
much  smaller  nondimensional  impulse  Zx,  in  this 
case,  as  shown  in  Table  1. 


The  predictions  of  the  analytical  model  are 
compared  with  explicit  finite  element  computations 
using  LS-DYNA.  In  the  computations,  the  cellular 
core  was  represented  by  a  single  row  of  solid 
elements  with  total  thickness  =5  cm,  using 
material  model  26  (*MAT_HONEYCOMB)  with 
Pq  =  250  kg/m^,  cTq  =  1  MPa,  and  =  0.7.  A 
large  elastic  modulus  of  E”  =  700  GPa  was  used  to 
represent  the  “rigid”  portions  of  the  idealized 
stress-strain  relationship  in  Fig.  1(b),  and  Poisson’s 
ratio  was  set  to  zero.  The  material  viscosity 
coefficient  p  was  set  to  0.001,  and  150  elements 
were  found  to  be  sufficient  for  convergence. 

The  front  and  back  faces  were  represented  in 
the  computations  by  added  nodal  masses,  and  two 
different  mass  distributions  were  considered,  as 
indicated  in  Table  1.  The  “pendulum”  case 
corresponds  to  the  blast  pendulum  experiments  of 
[1],  with  the  large  back-face  mass  representing  the 
pendulum.  The  “sandwich”  case  corresponds  to  the 
sandwich  plates  of  [3]  and  [5],  with  equal  front- 
face  and  back-face  masses. 

In  Figs.  2  and  3,  computational  results  are 
compared  with  predictions  of  the  analytical  model, 
and  good  agreement  is  observed.  Results  are 
plotted  against  nondimensional  time,  r  =  (cTq  /  )t . 
The  nondimensional  velocities  in  Figs.  2  and  3  are 


0  0.2  0.4  0.6  0.8  1 

Nondimensional  time,  x 


INFLUENCE  OF  MASS  DISTRIBUTION 

Fig.  3(a)  shows  contours  of  the  critical 
nondimensional  impulse  Zx  for  which  complete 
densification  of  the  core  is  first  achieved.  These 
contours  correspond  to  the  limiting  case  of  a  Dirac 
delta  impulse  (  Z^q  ^  oo  )  and  were  obtained  by 

numerical  solution  of  Eq.  (4).  Fig.  3(a)  shows  that 
increasing  the  mass  fraction  in  the  core  and  in  the 
front  face  increases  the  impulse  capacity  of  the 
sandwich  system.  However,  Fig.  3(b)  shows  that 
increasing  the  mass  fraction  in  the  core  and  in  the 
front  face  also  leads  to  increased  back-face 
accelerations,  thus  sacrificing  a  protective  function 
of  the  cellular  core.  The  nondimensional  back- face 
accelerations  presented  in  Fig.  3(b)  are  defined  as 
^2  -  {ml cTq )u2  .  It  follows  from  Eq.  (2)  that  the 
peak  back-face  accelerations  occur  at  the  instant  of 
complete  compaction  (x  =  0),  for  which 

^2  -  /  ^2  ^2  =  1  /  7/2 .  A  design  optimization 

problem  can  be  posed  by  seeking  to  maximize  the 
impulse  Zx,  that  can  be  sustained  while  limiting  the 
back-face  accelerations  to  a  specified  level.  Fig.  4 
shows  a  contour  plot  of  the  maximum  impulse  Zx, 
that  can  be  sustained  with  accelerations  limited  to 

-  5 .  The  grey  curve  in  Fig.  4  corresponds  to 

1/7/2  “  ^  •  Below  this  curve,  the  values  of 


Figure  2.  Comparison  of  LS-DYNA  computations  ( — )  with  predictions  of  analytical  model  (o):  Nondimensional  front- 
face  and  back-face  velocities  for  (a)  “pendulum”  case;  (b)  “sandwich”  case. 
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Figure  3.  Contours  with  varying  mass  distribution:  (a)  Critical  nondimensional  impulse  required  for  eomplete 
eompaetion  of  eore  (Pq  ^  (b)  Peak  nondimensional  baek-faee  aeeeleration  at  eomplete  eompaetion  of  eore. 


maximum  impulse  correspond  to  complete 
compaction  of  the  core  and  are  the  same  as  in  Fig. 
3(a).  Above  this  curve,  5  at  complete 
compaction,  so  only  partial  compaction  is 
permitted  and  the  values  of  maximum  impulse  are 
less  than  in  Fig.  3(a).  In  the  shaded  region  of  Fig. 
4,  defined  by  {rj^  +  772)”^  >  5 ,  a^>5  at  initiation  of 
compaction,  so  the  maximum  allowable  impulse  is 
zero.  It  is  evident  in  Fig.  4  that  for  a  given  mass 
fraction  in  the  core  ,  the  allowable  impulse  is 
maximized  along  the  grey  curve,  i.e.,  by  adjusting 
the  mass  distribution  so  that  the  acceleration  at 
complete  compaction  equals  the  allowable  value. 


0  0.2  0.4  0.6  0.8  1 

Fraction  of  face-sheet  mass  in  front  face 
m  /  ill  +  I2) 


Figure  4.  Contours  of  maximum  nondimensional 
impulse  with  nondimensional  hack-face 

accelerations  limited  to  d2  =  5. 
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Abstract.  In  this  study  we  examine  transient  stresses  in  elastic  and  viscoelastic  composite  strips. 
When  the  stress  is  the  result  of  a  stress  step  we  seek  combinations  of  material  parameters  that  lead  to  a 
minimum  peak  stress  in  the  layer  most  distant  from  the  applied  stress.  In  the  case  of  impact  by  a  rigid 
body  we  define  an  optimal  design  as  one  in  which  the  maximum  stress  is  achieved  in  the  layer  subject 
to  the  impact.  The  primary  purpose  of  these  results  is  benchmarking  larger  numerical  studies  coupling 
a  finite  element  code  (DYNA3D)  with  several  optimization  routines. 

Keywords:  Elastic  Waves,  Viscoelastic  Waves,  Optimization,  Inverse  Laplace  Transform. 

PACS:  62.30.+d,  43.20.Gp,  83.60.Bc 


INTRO 

Design  of  stress  wave  attenuators  is  a 
computationally  expensive  task  with  few  analytic 
solutions.  Therefore,  there  are  almost  no 
benchmarks  for  the  larger,  more  complex 
computations.  The  purpose  of  this  study  is  to  find 
reliable  solutions  to  transient  wave  propagation  in 
one-dimensional  composite  strips.  These  solutions 
are  examined  for  possible  optimal  designs. 

First,  we  introduce  the  idea  of  an  optimal 
design  by  considering  elastic  strips.  The  following 
sections  consider  viscoelastic  strips,  and  altering 
the  applied  stress  to  low-velocity  impact. 

ELASTIC  STRIPS 

For  comparison,  we  report  the  relevant  result 
from  a  homogenous  strip  of  an  elastic  material. 
This  will  be  used  to  benchmark  the  performance  of 
composites. 


Homogeneous  Elastic  Strip 

Consider  an  elastic  material 
occupying  0  <  v  <  Z ,  initially  at  rest,  and  fixed  on 
the  right  {x  =  L).  Upon  application  of  a  stress 
step  on  the  left,  stress  propagates  through  material; 
hereafter,  referred  to  as  a  strip.  The  mathematical 
model  of  this  is  the  Initial  Boundary  Value 
Problem  (IB VP): 

0  =  i/(x,0),  0  =  9^w(x,0)  ^2) 

aQH{t)  =  Ed  ^u{0j) 

0  =  u{Lj) 

where  u  is  the  displacement,  cJq  the  magnitude  of 
the  applied  stress,  H(t)  the  Heaviside  step 
function,  E  the  elastic  modulus,  and  p  the 
density  of  the  strip.  The  speed  of  stress  waves 
isc  =  -yjE/ p  . 
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Fig.  1  is  a  stress  history  at  two  points  along  the 
strip,  x  =  LI4  and3Z/4.  These  points  were 
chosen  because  they  will  be  the  midpoints  of  the 
layers  of  the  composite. 


Figure  1.  The  stress  history  measured  at  two  points  of 
the  homogeneous  elastie  strip. 

The  relevant  feature  of  Fig.  1  is  the  peak  stress. 
At  each  point  of  measurement  the  maximum  stress 


Two  Layered  Elastic  Strip 

Now  consider  a  strip  composed  of  two  layers. 
Layer  1  occupies  0  <  x  <  L  /  2  and  layer  2, 
L 1 2  <  X  <  L  .  All  quantities,  displacement,  stress, 
density  and  elastic  modulus,  will  have  their 
corresponding  layer  indicated  by  a  subscript.  The 
displacements  in  each  layer,  u.,  obeys  the  partial 

differential  equation  and  initial  conditions  in  (1). 
The  boundary  conditions  of  (1)  are  replaced  by 


We  restrict  our  attention  to  materials  where  the 
elastic  modulus  and  density  in  layer  2  are  scaled 
from  their  values  in  layer  1  by  the  same  factor,  a. 

E^=EJ  a,  p^=  pj  a.  (4) 

This  implies  =€2,  so  we  drop  the  subscript. 

Fig.  2  is  a  comparison  of  stress  histories  at 
x  =  3L/  4  for  the  homogeneous  strip  and  a 
composite  witha  =  1.1.  In  contrast  with  Fig.  1, 
the  stress  in  the  composite  exceeds  2(7q  . 


L 

Figure  2.  Stress  measured  at  x  =  3L/  4  for  (X  =  I 
(homogeneous  strip),  and  (X  =  1.1. 

This  phenomenon  is  common  for  almost  all 
composite  strips.  In  [I]  an  analytic  solution  is 
constructed  for  the  IB  VP  of  the  composite.  It  was 
shown  that  the  stress  amplitude  in  layer  2  will 
exceed  2(Jq  in  all  but  a  discrete  set  of  designs. 
These  ‘optimal’  designs  occur  when 


a^^H{t)  =  Ed  ^u^{xfi) 

0  =  U2{Lp) 


l  +  cos(;z'/^) 
l-cos(;z'/^)’ 


k  =  2,3,... 


(5) 


To  complete  the  model  we  need  interface 
conditions.  We  assume  that  the  displacement  and 
stress  are  continuous  at  the  layer  interface. 


i/j  {L  /  2-,  t)  =  U2{LI  2+,  t) 
E^d^u^  {L  /  2-,  t)  =  E2d^U2  {L  /  2+,  t) 


(3) 


For  optimal  designs  the  peak  stress  in  layer  2  will 
be2(jQ .  Fig.  3  shows  the  peak  stress,  in  layers  1 
and  2,  as  a  varies.  Optimal  designs  are  visible. 
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Figure  3.  Peak  stress  measured  at  two  points  of  the 
eomposite  strip  as  Ot  is  varied. 

In  the  following,  we  seek  results  similar  to  Fig. 
3.  We  will  consider  viscoelastic  strips  and  elastic 
strips  subject  to  low-velocity  impacts. 


VISCOELASTIC  STRIPS 

For  a  strip  composed  of  viscoelastic  layers  the 
IB  VP  (I)  must  be  adjusted.  We  keep  the  notation 
of  denoting  quantities  in  each  layer  with  subscripts. 
The  equation  of  motion  and  initial  conditions  are 

0  =  u.  (x,0),  0  =  3^1/ .  (x,0) 

The  boundary  and  interface  conditions  are 

0  =  ^2  (^5^)  ^7^ 

Wj [L / 2-, t)  2+, t) 

(Jj  {L  /  2-,  t)  =  a^{LI  2+,  t) 

The  complete  problem  statement  includes  the 
constitutive  relationship 

t 

a.  (x,  ^)  =  j  G.  {t  -  T)d^^u.  (x,  T)dT  ?  (8) 

0 


where  G-(t)  is  the  relaxation  modulus  of  layer/. 
We  restrict  our  attention  to  standard  linear  solids 

=  (9) 

G-  Q  is  the  initial  elastic  response  of  the  material, 
G-^  the  long  term  stress  response,  and  p.  the 
decay  parameter. 

We  make  the  following  parameter  choices. 
G2^q=EI(X 

^/,oo  “  /  ^5  Pi  ~  P\  ^  ^ 

The  speed  of  propagation  in  each 
layer,  c .  =  .^G.^q  /  p.  ,  is  the  same  and  will  be 

denoted  by  c . 

The  decay  parameter  is  a  viscosity  measure. 
Materials  with  small  P  will  behave  similarly  to 
elastic  materials  when  considering  transient 
stresses.  We  have  chosen the  same  in  each 

layer,  so  its  subscript  has  been  dropped. 

Equations  (6)  through  (9)  comprise  an  integro- 
differential  equation  for  the  propagation  of  stress 
waves  in  the  strip.  Due  to  its  complicated  nature, 
this  IB  VP  is  solved  numerically  using  the  Laplace 
transform  with  a  modified  DAC  algorithm  for 
inversion,  [2,3,4].  The  results  for  different  values 
of  a  and  three  values  P  are  displayed  in  Fig.  4. 


2.0  4.0  6.0 

cc 

Figure  4.  Peak  stress  at  x  =  3Z  /  4  as  61^  varies. 
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In  Fig.  4  we  see  that  as  a  varies  the  peak  stress 
in  layer  2  behaves  similar  to  the  elastic  composite. 
In  fact,  when  ^  is  small  the  results  are  almost 

indistinguishable  from  the  elastic  case.  As  P 

increases  the  viscous  nature  of  the  solid  becomes 
more  pronounced  and  the  difference  between 
optimal  and  nearby  designs  becomes  less  dramatic. 

LOW  VELOCITY  IMPACT 

Now  reconsider  the  two  layered  elastic  strip  in 
(1)  through  (4),  but  replace  the  stress  step  with  the 
impact  of  another  solid. 

Consider  a  rigid  body  of  mass  m  moving  to  the 
right  with  a  velocity  <c .  The  impact  is 
modeled  by  coupling  the  position  of  the  right  of  the 
rigid  ho&y, with  the  stress  in  the  strip 
through  the  condition 

(11) 

Fig.  5  contains  the  results  from  this  IB  VP  using  the 
characteristic  stress 

=Voa/^iA  •  (12) 

Optimal  designs  similar  to  those  in  Fig.  3  don’t 
exist.  However,  we  do  see  that  the  stress  in  layer  2 
is  the  lowest  stress  when  a  is  greater  than  1 .  This 
corresponds  to  a  composite  with  a  stiffer  layer  1 . 


a 

Figure  5.  Peak  stress  at  the  point  of  eontaet  and  the 
midpoint  of  eaeh  layer. 


CONCLUSIONS 

The  optimal  designs  of  composite  elastic  strips 
can  be  extended  to  strips  composed  of  viscoelastic 
materials,  with  similar  results.  When  low  velocity 
impacts  replace  stress  steps  the  optimizability  is 
lost.  However,  a  lower  peak  stress  is  recorded  in 
the  layer  2  when  the  front  layer  is  the  stiffer  of  the 
two.  Though  not  optimal,  this  is  still  a  desirable 
property  for  a  stress  wave  attenuator. 
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Abstract.  A  generalized  self-consistent  method  (GSCM)  is  used  to  evaluate  the  effective  moduli 
of  a  cracked  medium  and  the  wave  propagation  speed  within  this  medium.  Since  the  microcrack 
damage  will  in  general  have  an  angular  distribution,  the  wave  speed  will  be  anisotropic. 
Furthermore,  since  cracks  respond  differently  under  tensile  and  compressive  loads,  waves 
propagating  through  the  medium  under  different  loading  conditions  will  experience  different  wave 
speeds. 
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INTRODUCTION 

A  number  of  problems  of  physical  significance 
involve  waves  propagating  through  cracked  media. 
Examples  include  seismic  waves  traveling  through 
the  earth’s  crust,  the  impact  response  of  protective 
armor,  and  ultrasonic  damage  detection.  The 
microcracks  affect  both  the  static  properties  of  the 
medium  and  the  way  in  which  stress  waves 
propagate. 

When  the  wavelength  is  the  same  order  of 
magnitude  as  the  crack  length,  the  waves  will 
interact  with  individual  cracks.  There  are  many 
papers  in  the  literature  discussing  this  topic.  In  this 
study,  we  are  interested  in  the  problem  when  the 
wavelength  is  much  larger  than  the  crack  length, 
allowing  us  to  focus  on  how  the  damage  affects  the 
wave  speed  in  the  material. 

In  order  to  calculate  the  elastic  wave 
propagation  speed,  we  need  to  know  the  effective 
moduli  of  the  cracked  media.  Many  researchers 
have  proposed  different  methods  to  evaluate  the 
effective  moduli  of  the  cracked  media,  including 


the  Self  Consistent  Method  (Budiansky  and 
O’Connell  [1],  Horii  and  Nemat-Nasser  [2]),  Mori- 
Tanaka  [3]  method  (Benveniste  [4]),  Differential 
Scheme  Method  (Hashin  [5]),  generalized  self- 
consistent  method  (Aboudi  et  al.  [6],  Huang  et  al. 
[7],  Santare  et  al.  [8])  and  the  effective  self- 
consistent  method  (Zheng  et  al.[9]).  We  will  use 
the  anisotropic  GSCM  discussed  in  Santare  et  al.  [8] 

Sayers  and  Kachanov  [10,11]  studied 
microcrack  induced  elastic  wave  anisotropy.  These 
papers  connect  the  effective  moduli  of  a  cracked 
medium  to  the  elastic  wave  speed.  The  authors  use 
second-rank  and  fourth-rank  crack  density  tensors 
to  evaluate  the  effective  elastic  moduli.  Using  these, 
they  calculate  the  effective  moduli  for  perfectly 
aligned  and  randomly  oriented  crack  distributions. 
Following  these  results,  Schubnel  and  Gueguen 
[12,13]  used  micromechanical  and  statistical 
physics  to  evaluate  the  elastic  wave  velocities  and 
permeability  of  cracked  rocks. 

Most  of  the  above  results  are  only  applied  to 
randomly  oriented  or  perfectly  aligned  cracks.  And 
only  one  of  them  (Horii  and  Nemat-Nasser  [2]) 
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considers  the  difference  between  tension  and 
compression.  This  paper  will  investigate  the 
difference  between  elastic  wave  propagation  under 
tension  and  compression  loads  within  a 
microcrack-damaged  medium. 


Figure  1.  Four  load  cases 

ANISOTROPIC  EFFECTIVE  MODULI  OF 
MICROCRACKED  MEDIA 


region  of  area  A.  The  external  loads  pi,p2,h  and  q 
are  defined  in  Fig  1. 

Notice  that  if  the  loads  in  Fig.  1  are  made 
compressive,  most  cracks  will  not  open  in  Mode  I, 
thus  changing  the  energy  terms  in  Eqs.  (l)-(4).  A 
linear  solution  will  not  pick  up  the  physical 
difference  between  cracks  opening  in  tension  and 
closing  in  compression.  Consequently,  a  linear 
model  will  predict  the  same  moduli  in  tension  and 
compression.  To  account  for  crack  face  contact,  a 
finite  element  model,  with  nonlinear  contact 
elements  is  used.  This,  in  turn,  will  result  in  very 
different  effective  moduli  for  the  medium  when  it 
is  loaded  in  compression  than  when  in  tension.  For 
this  initial  study,  the  friction  associated  with  the 
contact  is  assumed  zero. 


Following  Santare  et  al.,  [8]  we  use  an 
anisotropic  Generalized  Self-Consistent  Method 
GSCM  to  evaluate  the  effective  moduli  of  the 
cracked  medium.  The  GSCM  equates  the  strain 
energy  in  a  damaged  region  to  the  sum  of  the  strain 
energy  in  the  region  if  there  were  no  damage  plus 
the  energy  associated  with  the  opening  of  the 
microcracks. 


In  a  general,  linear  2D  orthotropic  problem, 
there  are  four  independent  material  constants.  In 
order  to  evaluate  these  constants,  we  apply  four 
separate  load  cases  to  the  damaged  region  as 
shown  in  Fig.  1.  Applying  the  energy  balance  to 
each  of  these  cases  and  ignoring  the  effects  of 
friction,  results  in  four  coupled  equations  [8]  for 
the  corresponding  effective  moduli  , 

//*andir*. 


^  k=\  Q 

= ^pIa + 1;  J 


k=\C, 


1  1  M 

—q'-A  =  -q^A  +  ^] 


k=ic, 

M 


A  A  ;  1 

k=l  Q 


(1) 

(2) 

(3) 

(4) 


Here  E,  ju,  K  and  are  the  undamaged  moduli. 
\ui  ]  and  ^ .  are  the  load-specific  displacement  jump 


and  traction  along  each  of  the  M  cracks  in  the 


m 


-00  +00 

Figure  2.  Even  distribution  of  cracks 

In  both  tension  and  compression,  the  resulting 
effective  moduli  will  be  affected  by  two  properties 
of  the  microcrack  array,  its  crack  density  and  its 
orientation.  The  crack  density  is  characterized  by 
7]  =  Mc^/a  where  c  is  the  average  half  crack  length. 
The  orientation  is  characterized  by  a  probability 
distribution  around  a  mean  orientation.  Here,  the 
average  orientation  is  taken  to  be  parallel  to  the  x- 
axis  and  the  individual  crack  angle  is  0 .  The 
distribution  is  then  taken  to  be  constant  from 
to  /9q  as  shown  in  Fig.  2.  This  way  when  is  0 
degrees,  the  cracks  are  aligned  and  when  0^  is  90 
degrees,  the  cracks  are  randomly  oriented. 

WAVES  IN  DAMAGED  MEDIA 

Using  the  calculated  effective  moduli,  we  can 
predict  the  elastic  wave  speeds  in  cracked  media. 
In  order  to  find  the  wave  speed  of  plane  harmonic 
waves  in  direction  L  ,  we  need  to  solve  the 
following  eigenvalue  equation,  [14] 

det{4 -<y..pv"}  =  0-  (5) 
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(6) 


where,  A  =  L^CL  is  the  acoustic  tensor,  C  is  the 
elastic  tensor  and  L  is  a  3x2  direction  cosine 
matrix. 


[L]  = 


h  0 

0  L 


ih 


The  slowness  is  the  inverse  of  wave  speed: 


Figure  3.  Normalized  slowness  of  waves  in  eraeked  media  under  tension. 

The  arrows  point  to  the  direetion  where  0o  inereases.  (a)  ri=0.1,  (b)  ri=0.2,  (e)  ri=0.3,  (d)  ri=0.4 


Figure  4.  Normalized  slowness  of  waves  in  eraeked  media  under  eompression. 

The  arrows  point  to  the  direetion  where  0o  inereases.  (a)  ri=0.1,  (b)  ri=0.2,  (e)  ri=0.3,  (d)  ri=0.4 
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(7) 


ACKNOWLEDGEMENTS 


s<-'>(i)  =  ^ 

Fig.  3  shows  the  normalized  inverse  of  the 
wave  speed  (slowness)  in  the  cracked  media  under 
tensile  loads.  Fig.  4  shows  slowness  under 
compressive  loads.  The  outside  lines  are  for 
maximum  slowness  and  the  inner  lines  are  for 
minimum  slowness.  The  arrows  show  the  direction 
in  which  0o  increases  from  0  to  90  degrees. 

CONCLUSION  AND  DISCUSSION 

Microcrack  damage  has  a  significant  effect  on 
the  elastic  moduli  and  thus  will  change  the  wave 
propagation  within  the  media.  From  both  figures,  it 
is  evident  that  when  r|  increases,  the  wave 
slowness  increases  (speed  decreases).  This  is  due 
to  the  fact  that  as  r\  increases  there  are  more  cracks 
per  unit  area  and  the  medium  becomes  less  stiff  In 
addition,  when  the  angular  distribution  parameter 
00  decreases,  the  cracks  are  more  aligned. 
Therefore,  the  medium  becomes  more  anisotropic 
and  the  slowness  diagrams  become  less  circular. 

However,  we  can  see  that  when  0o=90°,  the 
figures  are  not  true  circles  which  means  the  results 
are  not  perfectly  isotropic.  The  reason  for  this  is 
that  in  the  analysis,  we  have  completely  decoupled 
the  tension  and  compression.  In  reality,  the  tension 
and  compression  moduli  have  an  effect  on  each 
other  through  shear  case  and  should  be  calculated 
simultaneously.  Currently,  we  are  working  on 
improving  the  analysis  to  address  this  issue. 

Comparing  Fig.  3  and  4,  we  can  see  that  as  the 
damage  increases,  the  wave  speeds  are  very 
different  under  tensile  and  compressive  loads.  As 
expected.  When  the  medium  is  loaded  in  tension, 
more  energy  is  associated  with  crack  opening 
(Mode  I).  Therefore,  the  modulus  decreases  more 
than  for  the  same  medium  loaded  in  compression. 
As  a  result,  under  tension  the  slowness  increases 
more  than  in  compression.  If  friction  were  included 
in  this  analysis,  it  would  affect  the  compressive 
results  more  than  the  tension  results  since 
compression  causes  more  crack  faces  contact. 
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Abstract.  The  cohesive  element  approach  is  getting  increasingly  popular  for  simulations  in  which  a 
large  amount  of  cracking  occurs.  Naturally,  a  robust  representation  of  fragmentation  mechanics  is 
contingent  to  an  accurate  description  of  dissipative  mechanisms  in  form  of  cracking  and  branching. 
This  paper  addresses  the  issue  of  energy  convergence  of  the  finite-element  solution  for  high-loading 
rate  fragmentation  problems.  These  results  provide  new  insight  for  choosing  mesh  sizes  and  size 
distributions  in  two  and  three-dimensional  fragmentation. 

Keywords:  dynamic  crack  propagation;  crack  branching;  fragmentation;  cohesive  element;  numerical 
convergence. 

PACS:  46.15.-x,  46.50.+a,  62.20.-x 


INTRODUCTION 

In  ceramic  materials,  it  is  known  that  strength, 
crack  initiation  and  propagation  are  affected  by  the 
presence  of  flaws  at  the  micrometer  scale.  Under 
dynamic  loading  conditions,  cracks  initiate  at  these 
flaws,  and  potentially  propagate  catastrophically  to 
cause  large-scale  structural  failure.  Multiple  crack 
initiate  at  seemingly  random  locations  and  material 
failure  occurs  through  a  complex  communication 
process  of  stress-wave  interactions  between  cracks. 
One  of  the  recent  numerical  approaches  for 
fragmentation  that  is  easy  to  implement  is  that  of 
cohesive  element  method.  Crack  branching  and 
fragmentation  are  natural  outcomes  of  the  method. 
In  the  creation  of  new  surfaces,  each  opening 
cohesive  element  dissipates  a  given  amount  of 
energy.  The  total  energy  dissipated  is  clearly 
related  to  the  crack  path  and  mesh  convergence.  It 
has  been  generally  observed  that  for  a  fixed  strain 
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rate,  finer  meshes  lead  to  higher  energy  dissipation 
due  to  more  microcracking  [1].  The  issue  of  energy 
convergence  during  the  cracking  process  is  of 
fundamental  importance  as  it  clearly  affects  the 
physical  robustness  of  the  numerical  methodology, 
and  is  the  main  purpose  of  the  paper. 

We  briefly  describe  the  finite  element 
methodology  including  the  cohesive  element 
approach  (illustrated  by  2D  crack-branching 
simulation)  and  analyze  convergence  issues  in  a 
simplified  one-dimensional  setting  and  the 
implications  for  extension  to  multidimensional 
finite  element  simulations. 

COHESIVE  ELEMENT  METHODOLOGY 

Our  adopted  framework  is  the  explicit  dynamics 
finite  element  analysis.  The  equations  of 
equilibrium  are  integrated  along  the  time  axis  by 
the  second-order  accurate  explicit  (Newmark) 
scheme.  For  stability  of  the  time  integration,  the 
time  step  At  must  be  lower  than  a  critical  value. 
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bistable,  which  is  related  to  the  dilatational  (the 
fastest)  wave  speed  and  the  (smallest)  mesh  size. 
We  use  the  initially  rigid  cohesive  law  [2]  and 
cohesive  elements  are  added  on  the  fly  (dynamic 
insertion).  The  initial  finite-element  mesh  is  free  of 
cohesive  elements,  and  as  the  dynamic  simulation 
proceeds,  cohesive  elements  are  inserted  at 
locations  where  the  stress  exceeds  a  critical  value 

CTc. 

In  addition  to  <Jc  (maximum  cohesive  force)  and 
5c  (critical  opening  displacement),  a  third 
parameter  Gc,  which  is  not  independent  of  the  other 
two  corresponds  to  the  area  under  the  <Jc  — 4  curve 
and  is  the  fracture  energy  that  is  needed  to  fully 
open  a  unit  area  of  crack  surface.  Resolution  of 
appropriate  length  and  time  scales  associated  with 
the  fragmentation  process  requires  proper  choice  of 
the  material  parameters.  Expressions  for  cohesive 
zone  length  a  have  been  obtained  by  [4,5].  The 
finite  element  mesh  size  must  be  smaller  than  a 
in  order  to  resolve  the  physics  of  dynamically 
failing  material  near  a  crack  tip.  It  has  also  been 
shown  [3]  that  cohesive  elements  introduce  a  time 
scale  into  the  problem  associated  with  the  opening 
of  an  isolated  micro  crack  which  determines 
whether  a  particular  externally  applied  strain  rate  is 
perceived  as  a  low  strain  rate  or  a  high  strain  rate 
during  the  dynamic  decohesion  process.  In  order  to 
properly  resolve  the  unloading  part  of  the  cohesive 
law  in  a  dynamic  fragmentation  simulation,  the 
integration  time  step  must  be  roughly  an  order  of 
magnitude  smaller  than  this  time  scale. 

Previously,  we  obtained  probing  crack 
branching  results  in  3D  PMMA  plates  [6].  For  the 
approach  in  2D,  similar  to  Falk  et  al  [7]  we  meshed 
a  h=3mm  square  block  of  PMMA  with  an  edge 
crack  of  0.25  mm.  Finite-element  meshes  of 
different  sizes  (smaller  than  the  cohesive  zone 
length  estimate  in  [5])  were  used. 

Crack  branching  results  were  obtained  for  all 
meshes  at  sufficiently  high  energy  levels  for 
different  loading  conditions,  indicating  the 
method’s  robustness.  However  an  increase  in  total 
cohesive  energy  dissipated  with  mesh  refinement, 
was  observed.  Closer  examination  of  the  cracked 
surfaces  showed  increased  microcracking  and 
crack  path  variation  with  mesh  refinement  (also 
observed  by  Ruiz  et  al  [1])  leading  us  to  question 


the  possibility  of  energy  convergence,  and  mesh- 
independent  microcracking.  These  questions  are 
addressed  in  the  next  section  for  an  idealized  ID 
setting. 

DYNAMIC  FRAGMENTATION  MODEL 


Figure  1:  Schematic  of  dynamic  ring  fragmenting  under 
explosive  loading 

For  simplicity  we  consider  the  most  elementary 
model  of  an  expanding  ring  (Fig.  1)  which  has 
been  discussed  in  earlier  work  [3,8].  We  assume  a 
linearly  elastic  material  and  insert  ID  initially-rigid 
cohesive  elements  upon  reaching  a  critical  stress  to 
monitor  multiple  cracks  nucleation,  growth,  and 
interactions.  The  ID  setting  allows  the  use  of  the 
method  of  characteristics  to  track  stress  wave 
interactions  in  the  elastic  material  and  permits  high 
mesh  densities  at  an  affordable  cost.  The  results 
were  also  compared  to  finite-element  calculations 
and  no  significant  differences  were  observed.  The 
ring  is  loaded  explosively  with  different  radial 
speeds.  Fragment  statistics  are  collected  when  all 
the  fragments  are  formed.  Fragmentation 
simulations  were  performed  for  different  mesh 
sizes  (from  128  to  1  million  nodes)  and  for  strain 
rates  ranging  from  5x10^  s'^  to  5x10^ 

Figure  2  represents  the  energy  dissipated 
(curves  in  solid  lines)  in  the  process  of  cohesive 
fracture.  This  includes  cohesive  elements  that  have 
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been  either  fully  or  partially  damaged.  We  observe 
a  non-monotonic  convergence. 


Figure  2:  Cohesive  energy  dependenee  on  mesh  size  for 
uniform  and  random  mesh  spaeing 

For  coarse  meshes  the  energy  dissipated  is 
proportional  to  the  number  of  nodes  in  the  mesh. 
This  strong  mesh  dependence  clearly  indicates  that 
these  mesh  sizes  under  resolve  the  fragmentation 
process  (all  the  cohesive  elements  are  fully  broken 
and  thus  the  fragment  size  is  mesh  size  dependent). 

The  cohesive  zone  size  estimate  is  shown  by 
the  vertical  line  in  Figure  2.  Clearly,  for  this  highly 
dynamic  problem,  this  estimate  of  mesh  size  is  not 
satisfying,  as  it  cannot  resolve  the  fragmentation 
process.  In  light  of  this,  one  may  question  the 
validity  of  using  cohesive  zone  length  scale 
estimates  in  any  dynamic  finite-element  problem. 
For  illustration,  consider  the  highest  strain  rate  case 
(5x10^  s'^).  We  note  that  the  dissipated  energy 
increases  up  to  lO"^  nodes.  Although  lO"^  nodes  is  a 
relatively  small  mesh  in  a  ID  setting,  its  2D 
equivalent  is  10^  nodes.  Many  published  results 
using  the  cohesive  element  approach  in  a 
multidimensional  space  have  been  obtained  with 
smaller  meshes.  We  postulate  that  observations  of 
energy  increase  with  mesh  refinement  are  due  to 
microcracking  being  not  properly  resolved, 
although  the  increase  may  not  be  linearly 


proportional  to  the  node  number  as  boundary 
conditions  may  deviate  from  the  simple  uniaxial 
state  of  stress  here. 

A  significant  finding  of  this  paper  is  that  for 
sufficiently  large  meshes  one  can  properly  resolve 
microcracking.  Smooth  convergence  is  observed 
beyond  the  peak  of  the  curve  This  is  encouraging 
as  it  indicates  that  although  cohesive  elements  are 
inserted  at  many  more  nodes  than  necessary  most 
are  not  severely  damaged  and  do  not  impact 
significantly  the  energy  balance  of  the  problem. 
However  it  is  disturbing  to  note  that  this  translates 
to  requiring  very  fine  meshes  (10^^  nodes  in  2D) 
beyond  any  available  computational  power. 

Using  a  smaller  time  step  did  not  affect  the 
non-monotonic  nature  of  the  curve  indicating  that 
perhaps  it  was  not  a  numerical  issue.  It  was  shown 
previously  [8]  that  a  simple  universal  law  (Weibull 
distribution)  captures  fragment  sizes  for  all  strain 
rates.  Nature  wishes  to  maximize  entropy  and  stays 
away  from  uniform  fragment  sizes  due  to  the 
inherent  randomness  associated  with  generating 
fragments  (which  uniform  meshes  do  not  capture 
well).  Uniform  meshes  severely  constrain  the 
fragmentation  event  and  its  desire  to  maximize 
entropy  and  seems  to  explain  the  observed  non¬ 
monotonic  convergence. 


Figure  3:  Mesh  dependence  of  fragment  size  distribution 
(random  mesh  spacing)  for  5x10^  s'^ 

The  curve  in  dotted  lines  in  Figure  2  represent 
the  results  obtained  for  random  meshes  obtained  by 
shifting  each  node  by  a  random  amount  in  between 
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±  0.4he,  where  he  is  the  mesh  size  in  the 
corresponding  uniform  mesh.  The  convergence  is 
now  monotonic  and,  remarkably,  it  is  up  to  two 
orders  of  magnitudes  faster.  In  addition  it  yields 
smooth  and  detailed  statistics  on  fragment  sizes 
(Fig.  3).  The  required  number  of  nodes  and  their 
equivalent  in  higher  dimensions  is  manageable 
with  parallel  simulations.  Numerical  modeling  of 
fragmentation  thus  has  an  odd  characteristic: 
uniform  meshes  should  be  avoided  at  all  cost. 


Instead  of  shifting  by  ±0.4he  we  tested  for 
varying  degrees  of  mesh  heterogeneity  and  found 
that  random  displacements  of  nodes  by  an  amount 
in  between  ±0.05he  suffice  to  improve  convergence 
by  roughly  two  orders  of  magnitude. 

As  discussed  earlier,  at  least  in  a  ID  setup, 
cohesive  zone  estimates  such  as  in  [5],  lead  to 
meshes  too  coarse  to  capture  all  microcracking 
events.  A  key  finding  of  our  previous  work  was 
that  a  universal  law  may  be  used  to  capture  the 
dependence  of  the  average  fragment  size  with 
strain  rate.  This  empirical  estimate  was  proposed  as 
an  alternative  to  Grady’s  energy-based  model  [9], 
which  does  not  take  into  account  wave  reflections. 
Using  the  average  fragment  size  relation  in  [10] 
and  Figure  2,  we  obtain  an  empirical  estimate  for 
required  mesh  size  (shown  as  dotted  line  in  Figure 
2)  as: 
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Using  this  estimate  and  perhaps  others 
incorporating  the  time  scale,  converged  2D 
fragmentation  simulations  are  within  reach.  3D 
simulations,  however,  remain  too  computationally 
intensive.  Rate  dependent  cohesive  laws,  for  which 
Gc  becomes  an  increasing  function  of  strain  rate, 
are  promising  in  this  regard  [11]. 

CONCLUSIONS 

Standard  cohesive  zone  size  estimates  under 
resolve  the  number  of  nodes  necessary  for  attaining 
mesh  independence  in  dynamic  fragmentation 
simulations.  Substantial  computational  evidence 
demonstrates,  for  the  first  time,  that  the  cohesive 


element  method  converges  in  an  energetic  sense. 
Microcracking  events  and  the  ensuing  fragment 
sizes  distributions  are  statistically  mesh 
independent  for  sufficiently  fine  meshes.  Random 
mesh  spacing  led  to  much  faster  and  smoother 
convergence  than  uniform  meshes,  and  that  even 
for  very  small  random  perturbations.  A  simple 
mesh  density  estimate  based  on  the  dependence  of 
the  average  fragment  size  on  the  strain  rate  and 
elastic  and  fracture  material  properties  is  proposed. 
This  estimate  provides  a  clear  roadmap  for  more 
complex  loading  conditions  including  biaxial 
loading  of  ceramic  plates. 
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3  SOUTHWEST  RSRCH  INST 
C  ANDERSON 

J  WALKER 
K  D  ANNEMANN 
PO  DRAWER  28510 
SAN  ANTONIO  TX  78284 

I  LIVERMORE  SOFTWARE 
TECH  CORP 
J  O  HALLQUIST 
2876  WAVERLY  WAY 
LIVERMORE  CA  94550 

I  TEXAS  A&M  UNIV 
DEPT  OF  MATH 
J  WALTON 

COLLEGE  STATION  TX  77843 


NO.  OF 

COPIES  ORGANIZATION 

5  UNIV  OF  NEBRASKA 
DEPT  OF  MECH  ENGRG 
D  ALLEN 
F  BOBARU 
Y  DZENIS 
G  GOGOS 
J  TURNER 
LINCOLN  NE  68588 

2  JOHNS  HOPKINS  UNIV 
DEPT  OF  MECH  ENGRG 
K  T  RAMESH 
FZHOU 

LATROBE  122 
BALTIMORE  MD  2I2I8 

I  WORCESTER  POLYTECHNIC  INST 

MATH  SCI 
K  LURIE 

WORCESTER  MA  01609 

4  UNIV  OF  UTAH 
DEPT  OF  MATH 
A  CHERKAEV 
E  CHERKAEV 
T  FOLIAS 
G  MILTON 

SALT  LAKE  CITY  UT  841 12 

I  PENN  STATE  UNIV 

DEPT  OF  ENGNG  SCI  &  MECHS 
F  COSTANZO 

UNIVERSITY  PARK  PA  168023 

3  UNIV  OF  DELAWARE 
DEPT  OF  MECH  ENGRG 
T  BUCHANAN 

T  W  CHOU 
M  SANT ARE 
126  SPENCER  LAB 
NEWARK  DE  1 97 1 6 

I  UNIV  OF  DELAWARE 

CTR  FOR  COMPOSITE  MTRLS 
J  GILLESPIE 
NEWARK  DE  1 97 1 6 

3  SRI  INTERNATIONAL 
D  CURRAN 
D  SHOCKEY 
R  KLOPP 

333  RAVENSWOOD  AVE 
MENLO  PARK  CA  94025 


NO.  OF 

COPIES  ORGANIZATION 

I  VIRGINIA  POLYTECHNIC  INST 
COLLEGE  OF  ENGRG 
RBATRA 

BLACKSBURG  VA  24061-0219 

I  COMPUTATIONAL  MECHANICS 
CONSULTANTS 
J  A  ZUKAS 
POBOX  II3I4 
BALTIMORE  MD  21239-0314 

I  KAMAN  SCIENCES  CORP 
D  L  JONES 

2560  HUNTINGTON  AVE  STL  200 
ALEXANDRIA  VA  22303 

I  APPLIED  RSRCH  ASSOC 
D  E  GRADY 

4300  SAN  MATEO  BLVD  NE 
STE  A220 

ALBUQUERQUE  NM  871 10 

I  INTL  RSRCH  ASSOC  INC 
D  L  ORPHAL 
4450  BLACK  AVE 
PLEASANTON  CA  94566 

I  AKT  MISSION  RSRCH  CORP 
M  EL  RAHEB 
23052  ALCALDE  DR 
LAGUNA  HILLS  CA  92653 

1  WASHINGTON  ST  UNIV 
SCHOOL  OF  MECH  &  MATE  ENGRG 
J  L  DING 

PULLMAN  WA  99164-2920 

2  WASHINGTON  ST  UNIV 
INST  OF  SHOCK  PHYSICS 
Y  M  GUPTA 

J  ASAY 

PULLMAN  WA  99164-2814 

I  ARIZONA  STATE  UNIV 
MECH  &  ARSPC  ENGRG 
D  KRAJCINOVIC 
TEMPE  AZ  85287-6106 


NO.  OF 

COPIES  ORGANIZATION 

I  UNIV  OF  DAYTON  RSRCH  INST 
NSBRAR 

300  COLLEGE  PARK 
MS  SPC  I9II 
DAYTON  OH  45469 

I  DIRECTOR 
USARL 
K  WILSON 
FRENCH  DEA  1396 
ADELPHI  MD  20783-1 197 

I  TEXAS  A&M  UNIV 

DEPT  OF  GEOPHYSICS 
T  GANGI 

COLLEGE  STATION  TX  77843 

I  UNIV  OF  SAN  DIEGO 

DEPT  OF  MATH  &  CMPTR  SCI 
A  VELO 

5998  ALCALA  PARK 
SAN  DIEGO  CA92II0 

I  NATL  INST  OF  STAND  &  TECH 

BLDG  AND  FIRE  RSRCH  LAB 
JMAIN 

100  BUREAU  DR  MS  861 1 
GAITHERSBURG  MD  20899-861 1 

I  BUCKNELL  UNIV 

DEPT  OF  MECH  ENGRG 
C  RANDOW 
LEWISBURG  PA  17837 

ABERDEEN  PROVING  GROUND 

73  DIR  USARL 

AMSRD  ARE  WM 
J  ALEHOUSE 
SKARNA 
J  MCCAULEY 
J  SMITH 
T  WRIGHT 
AMSRD  ARE  WM  B 
J  NEWILL 
M  ZOLTOSKI 
AMSRD  ARE  WM  BA 
AMSRD  ARE  WM  BC 
P  PLOSTINS 


NO.  OF 

COPIES  ORGANIZATION 


NO.  OF 

COPIES  ORGANIZATION 


AMSRD  ARE  WM  BD 
P  CONROY 
B  FORCH 
RLIEB 

R  PESCE  RODRIGUEZ 
BRICE 

AMSRD  ARE  WM  BF 
D  WILKERSON 
AMSRD  ARE  WM  EG 
E  SCHMIDT 
AMSRD  ARE  WM  M 
S  MCKNIGHT 
AMSRD  ARE  WM  MA 
R  JENSEN 

M  VANLANDINGHAM 
E  WETZEL 
AMSRD  ARE  WM  MB 
M  BERMAN 
L  BURTON 
T  BOGETTI 
W  CHOWDHURY 
W  DE  ROSSET 
W  DRYSDALE 
A  FRYDMAN 
D  HOPKINS 
L  KECSKES 
THLI 

M  MINNICINO 
B  POWERS 
J  TZENG 

AMSRD  ARE  WM  MC 
R  BOSSOLI 
S  CORNELISON 
M  MAHER 
W  SPURGEON 
AMSRD  ARE  WM  MD 
B  CHEESEMAN 
ECHIN 
B  DOOLEY 
C  FOUNTZOULAS 
G  GAZONAS 
J  LASALVIA 
P  PATEL 
J  SANDS 
CF  YEN 

AMSRD  ARE  WM  RP 
J  BORNSTEIN 
AMSRD  ARE  WM  SG 
T  ROSENBERGER 
AMSRD  ARE  WM  T 
AMSRD  ARE  WM  TA 
S  SCHOENFELD 


AMSRD  ARE  WM  TB 
P  BAKER 
R  SKAGGS 
J  STARKENBERG 
AMSRD  ARE  WM  TC 
R  COATES 
K  KIMSEY 
D  SCHEFFLER 
S  SCHRAML 
AMSRD  ARE  WM  TD 
S  BILYK 
T  BJERKE 
D  CASEM 
J  CLAYTON 
D  DANDEKAR 
M  GREENFIELD 
KIYER 
BLOVE 

M  RAFTENBERG 
E  RAPACKI 
M  SCHEIDLER 
S  SEGLETES 
T  WEERASOORIYA 
AMSRD  ARE  WM  TE 
J  POWELL 
B  RINGERS 
G  THOMSON 


NO.  OF 

COPIES  ORGANIZATION 


NO.  OF 

COPIES  ORGANIZATION 


I  DERA 

N  J  LYNCH 
WEAPONS  SYSTEMS 
BUILDING  A20 
DRA  FORT  HALSTEAD 
SEVENOAKS 
KENT  TN  I47BP 
UNITED  KINGDOM 

I  ERNST  MACH  INTITUT 
H  N AHAME 
ECKERSTRASSE  4 
D  7800  FREIBURG  I  BR  791  4 
GERMANY 

I  FOA2 

P  LUNDBERG 
S  14725  TUMBA 
SWEDEN 

I  PCS  GROUP 

CAVENDISH  LABORATORY 
W  G  PROUD 
MADINGLEY  RD 
CAMBRIDGE 
UNITED  KINGDOM 

I  CENTRE  D  ETUDES  DE  GRAMAT 
J  Y  TRANCHET 
46500  GRAMAT 
FRANCE 

I  MINISTERE  DE  LA  DEFENSE 

DR  G  BRAULT 
DGA  DSP  STTC 
4  RUE  DE  LA  PORTE  DISSY 
75015  PARIS 
FRANCE 

I  SPARE  DIRECTION  BP  1 9 

DR  E  WARINGHAM 
10  PLACE  GEORGES 
CLEMENCEOUX 
9221 1  SAINT  CLOUD  CEDEX 
FRANCE 

I  LMT  CACHAN 
J  F  MOLINARI 

61  AVE  DU  PRESIDENT  WILSON 
94235  CACHAN  CEDEX 
FRANCE 


I  TECHNICAL  UNIV  OF  CRETE 
G  EXADAKTYLOS 
DEPT  OF  MINERAL  RES  ENGNG 
CHANIA  CRETE 
GREECE 

I  ROYAL  MILITARY  COLLEGE  OF 
SCIENCE 

CRANEFIELD  UNIV 
J  MILLETT 

SHRIVENHAM  SWINDON 
SN6  SLA 

UNITED  KINGDOM 

I  UNIV  OF  MANCHESTER 
N  K  BOURNE 
PO  BOX  88 

SACKVILLE  STREET  MANCHESTER 
M60  IQD  UK 

1  BEN  GURIAN  UNIV  OF  NEGEV 
E  ZARETSKY 

DEPT  OF  MECH  ENG  BEER-SHEVA 
ISRAEL  84105 

2  RUSSIAN  ACADEMY  OF  SCIENCES 
INSTITUTE  FOR  HIGH  ENERGY 
DENSITIES 

G  I  KANEL 
S  V  RAZORENOV 
IVTAN  IZHORSKAYA  13/19 
MOSCOW  I274I2  RUSSIA 

I  INST  FUR 

MATERIALFORSCHUNG  II 
C  TSAKMAKIS 
POSTFACH  3640 
FORSCHUNGSZENTRUM 
KARLSRUHE 
D  76021  KARLSRUHE 
GERMANY 

I  TECHNICAL  UNIV  OF  DENMARK 

DEPT  OF  MATHEMATICS 
M  BENDSOE 
2800  LYNGBY 
DENMARK 

I  TECHNICAL  UNIV  OF  DENMARK 

DEPT  OF  MECH  ENGRG 
O  SIGMUND 
2800  LYNGBY 
DENMARK 


NO.  OF 

COPIES  ORGANIZATION 


I  NATL  TECH  UNIV  OF  ATHENS 
DEPT  OF  ENG  SCI 
I  VARDOULAKIS 
ATHENS  15773 
GREECE 

I  UNIV  OF  PATRAS 

DEPT  OF  CIVIL  ENGRG 
D  BESKOS 
PATRAS  26500 
GREECE 

I  ARISTOTLE  UNIV 

OF  THESSALONIKI 
DEPT  OF  MECH  &  MATES 
E  AIF ANTIS 
THESSALONIKI  54006 
GREECE 

I  DEMOCRITUS  UNIV  OF  THRACE 
DEPT  OF  CIVIL  ENGRG 
E  GDOUTOS 
XANTHI 
GREECE 


Intentionally  left  blank. 


